A method of generating initial conditions for cosmological N body simulations 



in 
o 
o 

(N 

c 
3 



M. Joyce, 1 D. Levesque, 2 and B. Marcos 2 

1 Laboratoire de Physique Nucleaire et de Hautes Energies, Universite de Paris VI, 
4, Place Jussieu, Tour 33 -RdC, 75252 Paris Cedex 05, France. 
2 Laboratoire de Physique Theorique, Universite de Paris XI, Batiment 210, 91405 Orsay, France. 

We investigate the possibility of generating initial conditions for cosmological N-body simulations 
by simulating a system whose correlations at thermal equilibrium approximate well those of cos- 
mological density perturbations. The system is an appropriately modified version of the standard 
"one component plasma" (OCP). We show first how a well-known semi-analytic method can be 
used to determine the potential required to produce the desired correlations, and then verify our 
results for some cosmological type spectra with simulations of the full molecular dynamics. The 
advantage of the method, compared to the standard one, is that it gives by construction an accurate 
representation of both the real and reciprocal space correlation properties of the theoretical model. 
Furthermore the distributions are also statistically homogeneous and isotropic. We discuss briefly 
the modifications needed to implement the method to produce configurations appropriate for large 
N-body simulations in cosmology, and also the generation of initial velocities in this context. 
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I. INTRODUCTION 

One of the major goals of current models in cosmology 
is to link the very small amplitude fluctuations observed 
in the cosmic microwave background to the large inhomo- 
geneities observed today in the distribution of galaxies. 
The former should serve as the initial seeds for the for- 
mation of the large structures observed, in the available 
time, through the attraction of gravity. In the absence 
of a full analytic understanding of the problem, gravita- 
tional N-body simulation has become an essential instru- 
ment in making this link (for reviews, see e.g. 0, 0, 
3] ) . In this paper we are concerned with one element of 
N-body simulation: the setting up of appropriate initial 
conditions (IC). 

The theoretical IC for such simulations are specified by 
the cosmological model considered, which gives the cor- 
relation properties, at the initial time of the simulation, 
of a continuous stochastic mass density field p{r). This 
field has a well defined mean density po(> 0), and fluctu- 
ations 8p{r) = p{r)/ po — 1, which are of small amplitude 
an generically assumed to be Gaussian. The correlation 
properties are thus fully specified by the power spectrum 
(PS) of the model. The subtlety in the generation of the 
IC for an N-body simulation is in the fact that a real- 
ization of this theoretical stochastic density field, which 
is a continuous function, must be discretised: it must 
be represented by an initial distribution in space of the 
N point particles used in the simulation. The correlation 
properties of these two qualitatively different kinds of dis- 
tributions, however, cannot be the same. At small scales, 
in particular, the discrete distribution has density fluc- 
tuations with infinite variance (points), while the contin- 
uous model has finite variance fluctuations at all scales. 
The correlation properties of the set of N discrete parti- 
cle must, nevertheless, approximate, in some appropriate 
sense, those of the theoretical model. 

There is one well known and universally employed 



method for generating the IC of N body simulations. It 
uses the so-called "Zeldovich approximation" 4, 5j. Par- 
ticles are displaced off a perfect lattice (or sometimes 
"glassy") configuration in a manner determined by the 
desired theoretical PS. Thus every realization of the the- 
oretical stochastic density field is mapped onto an initial 
configuration of the N particles. Each such configura- 
tion can thus be considered as a realizaion of a discrete 
stochastic density field. 

Usually the characterisation of the correlation prop- 
erties of this discrete distribution (i.e. of the generated 
discrete IC) are considered only in reciprocal space, i.e. 
in terms of a measured PS, which approximates well the 
PS of the input theoretical model for sufficiently small k 
(i.e. small compared to the Nyquist frequency character- 
istic of the discreteness of the underlying lattice). The 
question of the representation of real space correlation 
properties (e.g. mass variance in spheres and the two 
point correlation function) has been the subject of some 
discussion in the recent literature [27j ■ A detailed study 
of thi s qu estion is presented in a recent article by two 
of us [nj ■ The standard method is shown rigorously to 
give a distribution with a PS approximating well that of 
the input theoretical model below the Nyquist frequency 
of the lattice, for quite a broad class of input PS |28| . 
In real space, however, it can be shown that the corre- 
lation properties are generically a mixture, at all scales, 
between those of the underlying lattice and those of the 
theoretical model. The real space correlation function 
carries, notably, the trace of the underlying oscillating 
behaviour of the initial lattice correlation function. The 
mass variance in a sphere is also typically dominated by 
that of the lattice to scales which, depending on the in- 
put theoretical model and its normalisation, can be very 
much larger than the interparticle spacing. 

The question which inevitably follows is whether such 
a discrepancy between the correlation properties of the 
theoretical model and those of its discretisation are im- 
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portant for what concerns the interpretation of the re- 
sults of the dynamical evolution of the discrete N-body 
simulation under gravity. These evolved configurations 
are taken to represent, over a very wide range of scale 
(typically from well below the initial inter-particle sepa- 
ration up to scales a little smaller than the size of the pe- 
riodic box), the evolution of a continuous self-gravitating 
fluid with the initial theoretical correlation properties. Is 
it sufficient to have only a faithful representation of the 
initial correlation properties in reciprocal space? This is 
an open question, which we do not attempt to answer 
here. In order to answer it, however, it would evidently 
be useful to have other methods for generating discrete 
representations of the theoretical IC. More generally, for 
the study of gravitational evolution from near homoge- 
neous IC, it is of interest to have a diversity of methods 
for generating initial configurations. 

The aim of this work is to explore a possible such 
method which has been suggested in ^lj. The idea is 
simple: to find a system whose dynamics produces config- 
urations with the required correlation properties. Given 
that the IC describe very small fluctuations about ho- 
mogeneity, it is natural to seek a system whose thermal 
equilibrium gives rise to such configurations. In it 
has been noted that a simple modification of the inter- 
particle potential in the "one component plasma" (OCP), 
which is a well-known and much studied system in sta- 
tistical physics, gives rise to a PS at small wavenumbers 
which is identical to that of cosmological models, with 
fluctuations which are also Gaussian. Further it is argued 
qualitatively in [TT| that an appropriate modification of 
this potential at smaller scales could give rise to the cor- 
relations typical of standard "cold dark matter" (CDM) 
cosmologies (which are of the class currently favoured by 
cosmologists) . In this paper we show how a well known 
semi-analytic method, the hypernetted chain (HNC), for 
determining the correlation properties of such systems 
at equilibrium can be adapted to allow one to infer the 
required interaction potential given an input spectrum. 
We apply this inversion for some simple model CDM-type 
spectra, and find a potential with the qualitative features 
anticipated by [llj. We then verify the accuracy of this 
method with a simulation of the full molecular dynam- 
ics. This involves modifying the Ewald sum method for 
calculation of the potential. 

The paper is organised as follows. In the next section 
we discuss first the important underlying question of how 
the correlation properties of a continuous field may be 
represented by a discrete distribution. In section ITTT1 we 
briefly recall the basic features of the standard OCP, and 
explain the analytical statistical mechanics approach to 
the determination of its correlation properties. In par- 
ticular we describe the HNC approximation which gives 
a closed form for the determination of the two point cor- 
relation properties from the potential. We explain how 
it can be used in an inverted form to infer the required 
inter-particle potential for an input (physically reason- 
able) correlation function. We show then how this can 



be done for a CDM type spectra, giving numerical re- 
sults for some cases. In the next section we turn to the 
simulation of the molecular dynamics of these systems. 
We give results for several simulations which verify the 
results anticipated from the HNC, and provide actual 
realizations of discretisations of CDM type spectra. In 
the final section we discuss the possibility of using this 
method to generate configurations representing the IC of 
cosmological NBS, as in principle the method appears 
to work well and to provide an interesting, and perhaps 
more satisfactory, alternative to the standard algorithm. 
One issue is the generation of appropriate initial veloci- 
ties. We explain that this can be done as in the standard 
method, as the Zeldovich approximation gives a relation 
between the initial particle velocity and the gravitational 
field acting on it, which can also be applied to our con- 
figurations. A more non-trivial issue is that of the size 
of the configurations generated: the simulations we have 
performed to test the viability of the method are of a 
quite modest size — the largest simulation reported here 
has slightly less than 9000 particles — while current cos- 
mological simulations use, at the very least, hundreds of 
thousands of particles. Just as in the case of NBS itself, 
it will be necessary to speed up the algorithm developed 
here to attain this goal. The methods for this latter step 
already exist and can be applied to the more complicated 
two body potentials we study. 



II. REPRESENTATION OF CONTINUOUS 
SPECTRA WITH POINT DISTRIBUTIONS 

A. Discrete and continuous stochastic density fields 

Let us first recall (see e.g. E3) some basic prop- 
erties of the PS (or, in the terminology more common in 
statistical physics, structure function). We will use P(k) 
to denote this quantity when we refer to a continuous 
distribution, S (k) for the discrete case. For both cases it 
can be defined 29] as 

V^QG V 

where 6(k) is the Fourier Transform (hereafter FT) of 
the normalized fluctuation field (p(r) — po)/poj an d the 
angle brackets indicate an ensemble average. We will 
assume that the system is statistically isotropic so the 
PS will not depend on the orientation of the vector k 
i.e. P(k) = P(k). We will also assume [3(J statistical 
homogeneity (i.e. invariance of average quantities under 
translation). In this case the Fourier transform of the 
PS, for which we use the convention 
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is the reduced two point correlation function defined as 



B. Smoothing of discrete distributions 



m = Z(r) = 1 = <Wp(0)) ■ (3) 

The intrinsic difference between a continuous and dis- 
crete density field p{r) manifests itself in a qualitative 
difference between the mathematical properties of the 
two-point quantities in each case. In real space the cor- 
relation function £(r) has, for the class of finite one-point 
variance continuous fields which we consider, the prop- 
erty 



-1 < f(r) < £(0) < oo. 



(4) 



For the discrete case the one-point variance, which is 
equal to £(0), necessarily diverges because of the singular 
nature of the density field at any point. The correlation 
function can then be written 



£( r ) = J- 6 (r) 
n 



h(r) 



(5) 



where uq is the mean number density, S(r) is the (three 
dimensional) Dirac delta function, and h(r) is a non- 
singular function for all r which can be taken to have 
the property analogous to Eq.|@}. 

These properties in real space translate in k space into 
a difference in the asymptotic properties of the power- 
spectra at large k. The one-point variance of the density 
field is also given by the integral of the PS, and so for 
the continuous case we have 



lim k 3 P(k) = 



(6) 



in order that this variance be finite. In the discrete case, 
on the other hand, we have 



lim S{k) = — 

k— >oo 71q 

lim k 3 (S(k) ) =0. 

fe^oo no 



(7) 
(8) 



i.e. the divergence of the one-point variance is entirely 
associated to the "Poissonian" term in the PS, which is 
simply the FT of the delta-function singularity in real 
space explicit in Eq.JHJ. Note that both P(k) and S(k) 
are, by definition, positive functions, while S(k) — — 
is not. There is therefore no bound S(k) > 1/uq. In 
particular, one can have S(k) — > for k — > 0, in systems 
satisfying the constraint 



d 3 rh(r) 



1 

n 



(9) 



i.e. when there is appropriate anti-correlation to bal- 
ance the contribution to fluctuations at all scales from 
the Poissonian term associated to any discrete process. 
As discussed in , 0] , [l4| these correspond to highly 
ordered "super-homogeneous" systems. 



The intuitively evident fact that a discrete distribution 
can only represent the correlation properties of a contin- 
uous field above some scale — that characteristic of the 
"granularity" of the discrete distribution — is reflected 
mathematically in the differences just discussed between 
the properties in the two cases of the correlation function 
at small real space separations, and the PS at large wave- 
numbers. Let us suppose now that we have a discrete 
distribution with PS S(k), and a continuous distribution 
with PS P(k). What is meant when one says that the 
former is a discretisation of the latter? In what sense can 
we say that the former represents the correlation prop- 
erties of the continuous distribution with PS P(k)? The 
answer to this question is that there is in fact no unique 
prescription for passing between a discrete and continu- 
ous distribution. In particular taking formally the limit 
in which the number of particles goes to infinity at fixed 
mass density, which one might naively think to define 
the desired continuous limit, does not do so. Consider, 
for example, the case of an (uncorrelated) Poisson point 
process: as the number density is taken to infinity the 
fluctuations also go to zero. Thus the continuous limit is 
an exactly uniform distribution with P{k) = 0. 

As discussed in the most natural way of defining 
such a relationship is by an appropriate local smoothing 
i.e. we assume the represented density field is given by 
the convolution of the discrete distribution with a spatial 
window function W Rs (r) 



Pc(r)= / W Rs (\r-r'\)p d (f')d 3 



(10) 



where R s is the (single) characteristic smoothing scale 
and the realization of the discrete field is a sum over all 
the particles 



(11) 



and p c {f) is the corresponding realization of the contin- 
uous stochastic density field. We then have that 



P(k) = \W Ra (k)\ 2 S(k) 



(12) 



where W Rs (k) is the Fourier transform of ^(r). By 
the assumption that the window function gives a local 
smoothing, we mean that it is an integrable function. 
It is naturally normalised to unity (to conserve mass) so 
that Wr s (0) is equal to unity. Thus the PS of the discrete 
field must approximate well that of the continuous one 
for small k (i.e. k <C Rj 1 )- In real space the smoothing 
leads to the convolution relation 

€c(r)= / W Rs (r')W Rs (r")Ur + r' -r")d 3 r'd 3 r" 



(13) 

between the continuous correlation function £ c (r) and the 
discrete correlation function £d(r). One sees explicitly 
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how the singularity becomes regularised applying l|13|) to 
©: 



no 



-5(f) 



no 



W R Jr")W Rs (r + r")d 3 f' 



(14) 



Note that any pair consisting of a discrete and a con- 
tinuous density field, with PS P(k) and S(k) respec- 
tively, can be related to one another formally by Eq. (|12H , 
taken simply as a definition of the smoothing function 
[3l| . Whether S(k) can be considered to be a physically 
reasonable discretisation of P(k) depends then on the 
mathematical properties of this smoothing function i.e. 
whether it really represents a physical smoothing. It is 
useful, for what follows, to express the relation between 
the two spectra in a slightly different (but equivalent) 
form: 



S(k) = P{k) 



-D(k) 



(15) 



where no is the number density of the discrete distribu- 
tion, The function D(k) has then the properties imposed 
by Eqs.Q) and ©: 



lim D(k) = 1 

k—*co 

lim k 3 (D(k) - 1) = 0. 

k— >oo 

In real space one has analogously 

h(r) = Ur)-—FT[l-D(k)} 
n 



(16) 
(17) 

(18) 



where £ c (0 is the Fourier transform of P(k) i.e. the 
reduced two-point correlation function of the continuous 
model. Expressed in terms of the smoothing we have 
from Ea.lfnjl that 



\W Rs {k)\- 2 = l + 



n P(k) 



(19) 



Note that whether the smoothing which is associated to 
a D(k) is a physical smoothing depends, therefore, not 
only on its own properties, but also on those of P(k). 



C. Discretisation in the standard algorithm 

The standard method for generating IC for N-body 
simulations in cosmology superimposes correlated dis- 
placements on a "pre-initial" configuration of points rep- 
resenting "uniformity" . The latter is usually taken to be 
a perfect lattice, but sometimes also a "glassy" config- 
uration obtained by running the N-body code with the 
sign of gravity reversed (i.e. with Coulomb forces) and an 
appropriate damping term. Both of these configurations 
represent an unstable equilibrium point (approximately 
in the latter case) under their self-gravity, and evolve 
only when the perturbations are applied. The precise 



sense in which this method leads to a representation of 
the continuous model (and whether this is adequate for 
the dynamics which is simulated) is not discussed in the 
literature, but implicitly the criterion used appears to be 



S(k) « P(k) for k < k d 



(20) 



where kd is the wave-number characteristic of the dis- 
creteness (i.e. the Nyquist frequency in the case of the 
lattice). In ^(| we have investigated systematically the 
real and reciprocal space correlation properties of config- 
urations generated by this algorithm, using a formalism 
developed in [l^. We give, in particular, exact expres- 
sions for the PS S(k) of the generated distribution, for a 
given input theoretical PS P(k). From this one can infer 
explicitly the form of D(k), as defined in Eq. l|15|l . and 
the corresponding "smoothing" given through Ea. (|19fl . 
The essential point is that D(k) carries all the correlation 
of the "pre-initial" configuration. Indeed one has, to a 
first approximation, simply that D(k) « noPi n (k) where 
Pin(k) is the PS of the "pre-initial" (lattice or glass) 
configuration. In fc-space this correlation is "localised" at 
large k, so that Ea. (|2U[l is satisfied. In real space, how- 
ever, the second term on the right hand side of Eq. i|18l 
is not localised to below the inter-particle distance, and 
completetly dominates the small amplitude correlations 
of the input model |33- Correspondingly the algorithm 
generically does not give a discretisation which can be 
interpreted as a localised smoothing in the sense we 
have discussed above: the function Wr s (r) determined 
through Ea. lfT^I (or, equivalently, Ea. l(T§|) is also highly 
delocalised). This is simply because the underlying reg- 
ularity of the "pre-initial" configuration, which is not a 
feature of the continuous model, cannot be removed by 
a localised smoothing. 



D. Determination of the PS of a new discretisation 

We investigate here a different method for discretising 
a given input PS P(k). The principle is to seek to gen- 
erate a distribution with an S(k) given through Eci. (|15fl . 
where for D(k) we will choose a smooth function of k, 
characterised by a single scale kd, and interpolating be- 
tween zero for k < kd and unity for k > kd (and in keep- 
ing with the asymptotic properties required Eas. i|ltj|) and 
lEl). The scale kd will be chosen of order the inverse of 
the mean particle separation a (see below for the exact 
definition we use). Further the function D(k) will be such 
that the FT of (D(k) - 1) in Eq.lO is localised strongly 
in real space on the scale a. Thus, by construction, we 
will converge to 



S(k) sa P(k) for fc<fc rf 
h(r) « £c( r ) for r a 



(21) 
(22) 



As we have noted, whether this choice of D(k) corre- 
sponds to a physical smoothing, in the sense we have 
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discussed, depends also on the properties of P(k). For 
the well-behaved P(k) we will consider we expect this to 
be the case, but we will check explicitly that the function 
Wr s (r) is smooth and integrable. 

The precise scale k < kd at which Eq. I|21|l holds will de- 
pend both on D(k) and on the form and normalisation of 
the PS. In the cosmological context P(k) is generically a 
monotonically decreasing function over a wide range of k 
for k < kd, and thus the dimensionless quantity noP(kd) 
gives a parametrisation of the relative amplitudes of the 
"continuous" and "discrete" parts of the full PS S(k). In 
the simulations of molecular dynamics described below 
we will take noP(kd) ~ 1. Thus we will have in this case 
Ea.(|2TTl for all k<k d , and (we will verify) Eq. lj22|) from 

In our explicit examples of the construction of S(k) we 
will make the simple choice D(k) = 1 — e~ fe / 2fe d, which 
evidently has the required asymptotic properties. It is 
important to note that we have not shown that the S(k) 
then given by Eq. I|15fl and such a choice of D(k) is neces- 
sarily the PS of a real discrete distribution |2J|. Indeed 
it is easy to see that the ansatz for S(k) may be unreal- 
izable in a discrete distribution: we have noted that the 
two-point correlation function h(r) of the discrete dis- 
tribution must satisfy by definition h(r) > — 1. Taking 
Ea. (|18|l . it is not difficult to verify that this condition 
places an upper bound on kd, of order the inverse of the 
average inter-particle distance |35j . Physically it is very 
reasonable that such a bound arises: taking kd larger 
than the inverse of the inter-particle separation one is re- 
quiring the discrete distribution to mimic the correlation 
properties of the continuous model in a regime where the 
intrinsic difference in the nature of the distributions is 
important. 



III. SEMI-ANALYTIC DETERMINATION OF 
THE POTENTIAL: THE HNC 

In this section we discuss first the OCP and a stan- 
dard method use to determine its correlation properties. 
The inversion of this method to determine the potential 
which would be expected to give reasonable desired input 
correlation properties is then explained. 



A. The standard OCP 

The OCP (for a review, see ^(J) is a system of posi- 
tive charged point particles ( "ions" ) interacting through 
a Coulomb (i.e. repulsive I/r) potential, and embedded 
in a uniform (rigid, non-dynamical) negatively charged 
background [36| . The latter gives overall charge neutral- 
ity, and a high degree of stability to the system. The sys- 
tem exhibits two phases at thermal equilibrium, a fluid 
phase and a solid phase. We will treat it always at den- 
sities and temperatures where it is in the fluid phase. In 
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FIG. 1: Correlation function of the OCP with interacting 
Coulomb potential for different temperatures (recall that F ~ 
1/T). 

this range of densities and temperature it can be consid- 
ered as completely classical. 

The equilibrium thermodynamics of the OCP is deter- 
mined by a single parameter, and not by its temperature 
and density independently. Because of the scale-free na- 
ture of the power-law interaction potential, there are only 
two characteristic length scales. One is specified by the 
number density, and is conventionally taken to be the 
"ion-sphere" (or simply ionic) radius a defined by 

/ 3 \ 1/3 

-(55) (23) 

where n = N/V is the number density of the N points in 
a volume V. The other scale is given by the distance at 
which the potential is of order the mean thermal kinetic 
energy. It is the dimensionless ratio of these two scales 
which parametrises the one dimensional phase space of 
the system at thermal equilibrium. Conventionally this 
parameter is taken to be 

T = P(Ze) 2 /a. (24) 

where (3 = l/(fcsT) and Ze is the ionic charge. It is ref- 
ered to as the "plasma parameter" (or simply "coupling 
constant" ) . 

Numerical calculations are necessary to compute reli- 
ably the correlation properties of the OCP. The correla- 
tion function g(r) (defined as g(r) — h(r) + 1) and PS 
S(k) are shown in Figs. ^ and [21 for different values of 
the coupling T. 

In the weak coupling limit (T < 1 i.e. high temper- 
ature/low density) the OCP can be shown to have two 
point correlations given approximately as 

Mr) = (g(r) - 1) (3- , k 2 = 47rZ 2 e 2 /3n. (25) 

r 

i.e. exponentially decaying anti-correlations. This is sim- 
ply a manifestation of screening, and k^ 1 is the Debye 
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FIG. 2: Power spectrum of the OCP with interacting 
Coulomb potential. Note that this is a linear-linear plot. 



length: 



\ D = (4TrZ 2 e 2 nf3) 



-1/2 _ 



(26) 



Particles are repelled locally around a given particle in 
such a way that the effective potential due to a particle 
decays exponentially in this way. At larger values of T 
(i.e. lower temperature/higher density) one sees that the 
correlation function develops a "bump" at small scales, 
indicating that the first neighbour is becoming increas- 
ingly localised. As T increases further several "bumps" 
develop (corresponding to first, second, third neighbours) 
which give to the correlation function and PS an oscil- 
latory structure, fore-shadowing the transition to the or- 
dered solid phase at T w 180 (for more details, see 0). 

By taking the Fourier transform it is easy to verify 
that, at small k the PS has the form 



(27) 



The behaviour S(k — > 0) — * we see in all cases is a 
manifestation in k space of the very strong suppression 
of the fluctuations at larger scales in the system due to 
screening. The fluctuations in particle number (i.e. in 
charge) in a volume are proportional to the surface of 
the volume, which is in fact the limiting behaviour pos- 
sible for any point distribution This is actually 
0, 0] a property shared by the PS of the canonical 
model in cosmology for "primordial" fluctuations (the so- 
called "Harrison- Zeldovich" (HZ) spectrum), which has 
the different small k behaviour P{k) ~ k. 

In it has been pointed out that, for a modification 
of the OCP with a repulsive potential v(r), screening 
generically leads to the behaviour for small k 



S(k) ~ FT[«(r)], 



(28) 



(where the FT is taken in the sense of a distribution). 
While for the Coulomb potential one obtains the be- 
haviour of Ea. H27|) . it is evident that by taking instead 



an interaction potential v(r) ~ 1/r 2 , one should obtain 



5(fc) ~27W f ° rfc<< £- 



(29) 



As we will discuss in further detail at the appropriate 
point below, cosmological PS which are used as the in- 
put for numerical simulations of structure formation are 
more complicated than the simple HZ form. In principle, 
however, it seems possible that an appropriate further 
modification of the interaction potential could allow a 
system at thermal equilibrium to produce more compli- 
cated PS. 



B. Determination of h(r) using the HNC equation 

The Hypernetted Chain Equation, which can be de- 
rived in a perturbative treatment of the full equilibrium 
thermodynamics (see e.g.^j|, chapter 5), gives a simple 
relation between the full two point correlation function 
and the interaction potential: 



h(r) = exp[— {3v(r) + h(r) — c(r)] — 1 . 



(30) 



The function c(r) is the direct correlation function. It 
gives the contribution to the correlation coming only 
from the interaction of two particles at the given sep- 
aration (i.e. neglecting all indirect correlation due to 
interaction with other particles). It is related to the full 
correlation function through the Ornstein-Zernike (OZ) 
relation 



h(r)=c(r)+n / dr'c(\r - r'\)h(r '). 



(31) 



Given the potential v(r), Eas- H^Ufl and i|31|) give a closed 
set of equations for the correlation function h(r), which 
can be solved by iteration as follows. Transformed to 
Fourier space Ea. l|31|) is a simple convolution product 

h(k) = c(k)+nh(k)c(k). (32) 

which allows one to express h{k) in terms of c(k) as 



h(k) 



1 — nc(k) 



It is convenient to define j(r) = h(r) 
FT j(k) is given in terms of c(k) as 



(33) 

c(r), of which the 



1 — nc{k) 



c(k). 



(34) 



We start with a first guess for c(r), denoted Co(r). One 
can take Co(r) = Vr, or, for example, co(r) ~ —(3v(r), 
which is the solution to the HNC equation Ea, 13UI) at 
leading order in expansion of the exponential. We can 
then use a Fast Fourier transform (FFT) to calculate 
co(fc), which then gives 70 (k) through l|34|) . With an in- 
verse FFT we find 70 (r), and then use the HNC equation 
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Eq. (|3"0|l to compute Ci (r) (using 70 (r) in the exponent to 
obtain Ci(r) + 70 (r) on the left hand side). The iteration 
process then proceeds with the computation of 7j(fe) with 
(|34[) . To ensure convergence, successive approximations 
on 7(7*) need to be taken, so the ith input is mixed lin- 
early with the precedent one: 



7-(r) = aji-i(r) + (1 - a)ji(r) 



(35) 



where < a < 1. The new 7^(r) is substituted in equa- 
tion Q34Jl to get Ci+i(r) and so on. In the simulations 
described below we took a — 0.5 which gives rapid con- 
vergence (less than one hundred iterations were neces- 
sary in all cases) . If there are problems with convergence 
(which can occur e.g. at larger densities) a value of a 
closer to 1 is taken. 

There is one additional elaboration of this method 
which is necessary when the potential is long-range, as 
it is for the case of the standard OCP and all the mod- 
ifications which we will consider in what follows. Since 



c(k) 



h{k) 



l + nh(k) n 



1 - 



n a S(k) 



(36) 



we have that c(k) diverges for k — ■+ 0, which is problem- 
atic numerically. This is dealt with in an analogous man- 
ner to that described in the Sect. II VI below for the calcu- 
lation of the force by the Ewald sum. One breaks c(r) into 
the sum of a short-range part c s (r) , with an analytic FT 
at k = 0, and a long part f(r), which contains the diver- 
gence in the FT. A typical long range part is u(r)erf(?7r) 
or v(r)(l — exp(— r/r)), where rj is a free positive parame- 
ter (on which the final result does not depend). The to- 
tal correlation function h(r) has no divergence, and thus 
7(r) is divided in the same way, j(r) — 7«(r) + /(r), with 
7 S (r) = h(r) — c s (r). The potential likewise is separated 
into a short and long range part j3v s {r) = v(r) + f(r), so 
that the HNC reads 



h(r) = exp[-/3v s (r) + 7«(r)] - 1. 



(37) 



When we compute Ea. H34(l we use the FT of the long- 
range part of the potential: 



c s {k) + f{k) 
l-n(c.(fc) +/(*)) 



c s (k). 



(38) 



All the computations are then done as described above 
but with c s (r) and 7 s (r) instead of c(r) and 7(7*), and 
using Eq. (J^J instead of Ea.p^l. 



Starting from an input model specified by a given PS 
S(k) we need just to calculate h(r) and c(r) (using the 
OZ relation Eq.J32J). This can most conveniently be 
done using FFTs. 

As noted above, when we treat the case of a PS with 
S(k — > 0) = 0, characteristic of a long-range interaction 
potential, we have a divergence at k = in c(k). Just as 
in the direct use of the HNC we deal with this numerically 
by dividing c(k) into two parts. The short-range part, 
which is regular at k = 0, can be taken to be 



^ = — ( 1 — hi* 

n \ n Q b(k) 



erfc(fc?y) \ 
n S {k) J 



(40) 



where Sa(k) is the functional form of S(k) at small k, and 
as above, 77 is a parameter on which the final result does 
not depend. The subtracted divergent piece is chosen (if 
possible) so that it can be Fourier transformed analyt- 
ically, and the full potential can thus be reconstructed 
easily from a determination of the short-range part of 
the potential from Ea. (|39|l using c s (r): 



(3v(r) — (3v s (r) — FT[c;(fc)], 



(41) 



Inversion of the HNC 



where 5;(fc) is the long-range part of c(fc), which corre- 
sponds to f(r) in the precedent subsection. 



IV. DETERMINATION OF THERMAL 
EQUILIBRIUM PROPERTIES USING 
MOLECULAR DYNAMICS 

The two numerical methods used widely in statisti- 
cal physics to study systems at thermal equilibrium are 
molecular dynamics and Monte Carlo simulations. We 
will discuss the former method, in which one evolves nu- 
merically the 3iV classical coupled equations of motions 
of a system of N interacting particles in a volume V. 
Finite-size effect are treated using periodic-type bound- 
ary conditions. We discuss the modification of a standard 
OCP code, developed by one of us (DL), to study the 
reliability of the (approximate) HNC method described 
in the previous section. We are interested therefore in 
simulating in this way a system of particles interacting 
through two-body potentials expected to generate sys- 
tems with two-point correlation properties like those of 
cosmological models. Further, if this kind of method is to 
be used to generate initial conditions for N body simula- 
tions of gravity, these simulations provide an instrument 
for producing the required configurations. 



It is simple to use the HNC equation in the inverse 
direction i.e. to determine an interaction potential v(r) 
which should give at thermal equilibrium desired two- 
point correlation properties, for we have directly from 
Eq.ljlUl that 

f3v(r) = h(r) - c{r) - \n[h(r) + 1]. (39) 



A. Discretisation of the Newton equations 

To discretise the equations of motion we use the Verlet 
algorithm. Performing a Taylor expansion of the position 
of a particle at times t + At and t — At about its position 
at time t, the position of the i-th particle is given to order 
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C((Ai) 4 ) by: 



fi(t + At) = 2f(t) - r(t - At) + 



(At)' 



N 



This algorithm, which is historically one of the earliest 
ones, has three important properties: it conserves energy 
very well, it is reversible (as the Newton equations) , and 
it is simplectic (i.e. it conserves the phase space volume). 
More refined algorithms have been proposed and used, 
but they often have less good conservation of energy at 
large times. Furthermore, the rapidity of the execution 
of the program is not determined by the computation of 
the new positions but by the calculation of the forces. 



B. Force calculation using the Ewald sum 

In our simulations N particles are placed in a cubic 
box of size L. To compute the interaction between the 
particles we apply the image method to minimize bor- 
der effects: an infinite number of copies of the system 
is created and the potential is computed considering not 
only the particles situated in the original box but also the 
particles of all the copies. Then if the particle i has co- 
ordinate fi, its copies will have cordinates fi + ftL, where 
ft is a vector with integer components. For a power-law 
interaction potential v(r) = r~ a the potential is then 



-nL 



(43) 



:)•» 



where qj is the charge of the particles and the asterisk de- 
notes that the sum ft — does not include the term i — j. 
In a numerical calculation the infinite sum Ea . I|43|) must 
be truncated. For a > 3 the potential is short-range and 
the approximation to compute the interaction potential 
between the i and j particles by taking only the interac- 
tion between i and the closest image of j is very good. 
When the potential is long-range (a < 3) this approxi- 
mation is no longer good, and indeed the sum appears 
to be formally divergent. For the case of the Coulomb 
potential, the presence of the neutralising uniform back- 
ground ensures that the potential of the infinite periodic 
system is well defined. A natural way of writing the sum 
in an explicitly convergent way taking this regularisa- 
tion into account is to separate the potential into a short 
range and long range part by introducing a parameter- 
dependent damping function f(f;a): 



Qj 



m. 



nL: a) 1 
+ — 



fin 



ftL; a) 



ftL\ a 



■ nL\ a 



(44) 

The first term on the r.h.s of Eq. 144fl is short-range and 
the second term is long-range. The procedure used in the 
Ewald summation method is to compute the first term 
in real space and the second in Fourier space. If the 



parameter a is appropriately chosen the real part con- 
verges well taking only the sum over the closest image, 
and the part of the sum in Fourier part is rapidly conver- 
gent. Physically the first term corresponds to a smearing 
of the original distribution, and the second term to the 
original point distribution sorrounded by a countercharge 
smeared distribution. Of course the sum of the two terms 
yields the original particle distribution. We write the po- 
tential energy then as: 



(45) 



Further it is convenient to separate out the zero mode in 
the long range part, writing 





k 




k=0 



. 

k^O 



(46) 



The function /(r; a) is chosen in the Ewald summation so 
that <p~ and & r e both rapidly convergent, and with 
a known analytical expression for its Fourier transform. 
The value of the term k = depends on how precisely the 
infinite sum in Eq. (|43|l is defined, and, as we will see fur- 
ther in particular examples, it is equal to zero in the pres- 
ence of the background because of the charge neutrality. 
This method of evaluating the potential energy using the 
Ewald Sum has been generalised for gen eric r~ a poten- 
tials [2(J , and for a Yukawa potential |2lj . In principle it 
may be used for other potentials. Note in particular that 
the Ewald method is applied to sum the long-range part 
of the potential: it remains valid if one introduces any 
additional short-range potential which can be absorbed 



without modification of 



/,(') 



We now give more 



detail first on its implementation for the standard OCP, 
and then for the potentials which will interest us here. 



1. Case 1: Coulomb interaction 

The /(r; a) function is usually chosen to be 

f(r: a) = erfc(a|f + ftL\) (47) 

where erfc is the complementary error function, erfc(x) 
1 — dt exp(— t 2 ). It is equivalent to smearing the 

charge distribution to obtain 

JV 
3=1 ft 

and introducing in Fourier space the original distribution 
plus the opposite smeared distribution. With this choice 
the short-range interaction energy is given by 



j=l n 



Ifji + HL\ 
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and the long-range part by 

JV 



(*) = §EE*p~ f ~" 2 



3=1 k^O 



F 6XP V4a 2 



cos(kfij). (50) 



The k = term is thus well defined only when the con- 
straint 



1 N 



(51) 



is satisfied i.e. only for a neutral distribution of charge. 
In the present case this is imposed by the uniform rigid 
negatively charged background assumed. 

An appropriate choice of a is a w 5.6/ L, where L is 
the size of the box. This gives good convergence in both 
(I49|l and l|5UI) . i.e. it includes only the first term n = 
in the first equation and not too many k in the second. 



the long-range part that comes from the subtracted di- 
vergence on the r.h.s. of Ea. H4U[) . the long-range part is 
thus in this case 



1.(0 



i =1 fc#o 



(55) 

where So(k) = Nk gives the small k behaviour of S(k). 
The real part of the potential is then: 



2ir 2 n 2 Nj3r 2 ' 



(56) 



Note that the parameter a in the Ewald sum needs to 
have the same numerical value as the parameter 77 in the 
HNC. 



V. GENERATION OF DISCRETISATIONS OF 
COSMOLOGICAL SPECTRA 



2. Case 2: l/r 2 potential 

With a potential l/r 2 a convenient choice for f(f;a) 
is M: 



/(r; a) = exp(a 2 |r + nL\ 2 ). 
The short-range part of the energy is 



EE 



iV 



9j- 



3=1 n 

and the long-range part 

2 JV 



exp(a 2 |Fy +nL\ 2 ) 
\r,i +nL\ 2 



(52) 



(53) 



= 7T EE *4 crfc ( ^ ) cos (^)- ( 54 ) 



3=1 k^O 



We will use the same value for a as in the Coulomb case. 
With this value of a the real part still converges rapidly 
and the Fourier part is much more rapidly convergent. 



3. Case3: inversion of HNC 

In what follows we will wish to simulate the molecu- 
lar dynamics of particles interacting through the poten- 
tial determined by the inversion of HNC equation as de- 
scribed in the previous section. As discussed in Sect. [HJ 
the small k behaviour of cosmological PS (the HZ spec- 
trum of perturbations), requires a long-range l/r 2 poten- 
tial. In the determination of the full potential through 
the inversion of the HNC, this piece is separated out by 
construction and the result is written as a sum of it and 
the short-range part subsequently determined. Taking 



N-body simulations of the formation of structure in the 
distribution of matter at large scales start from an initial 
time which is "recent" in terms of cosmological history. 
The universe has entered the phase in which its energy 
density is dominated by massive particles, and the evolu- 
tion of perturbations in the distribution of these particles 
at the scales considered is well approximated by Newto- 
nian gravity. The fluctuations at this initial time are still 
of small amplitude at the relevant physical scales, and 
the simulation follows this evolution through to today 
when very high amplitude fluctuations have formed at 
scales comparable to those on which they are observed 
to exist today. These initial conditions for simulations 
are generically Gaussian in current cosmological models, 
and thus fully specified by their PS. This PS is the result 
of the evolution up to this time, which can be calculated 
precisely in a given model (and depends on the various 
parameters characterising it) of the "primordial" fluctu- 
ations, which usually have the form given by so-called 
"scale- invariant" fluctuations [33- Because the fluctua- 
tions evolve in a non-trivial way for a finite time (un- 
til the time of "equality" , after which matter dominates 
over radiation) the resultant PS corresonds to the "pri- 
mordial" spectrum P(k) ~ k only up to a characteristic 
wave-number kt, above which it "turns over" to a differ- 
ent behaviour, with a PS which decreases as a function of 
k but with a functional behaviour which depends on the 
model. We will consider here the class of "cold dark mat- 
ter" (CDM) models which are those currently favoured 
as viable models to explain the diverse observations of 
large scale structure. In N body simulations the initial 
spectrum is usually taken in this case as given by a simple 
functional fit to the results of a numerical determination 
of the PS (see e.g. @): 



P(k) 



jV{z)k 



[1 + (aq + {bqy/ 2 + {cq) 2 



.2/v 



(57) 
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where q = fc/A is a rescaling of k by a dimensionless pa- 
rameter A which depends on the parameters of the CDM 
model (A = 0.21 for "standard" CDM), and v = 1.13. In 
units of hr 1 Mpc, where h is the Hubble constant today 
in units of 100 km/s/Mpc, one has a = 6.4, 6 = 3 and 
c = 1.7. The factor N(z) gives the overall normalisation 
of the PS, which is a function of the initial red-shift z (for 
a red-shift chosen in the matter dominated era, during 
which the fluctuations are, to a very good approxima- 
tion, simply amplified in the same way at all scales.) It 
is in principle fixed by the amplitude of fluctuations mea- 
sured in the cosmic micro- wave background (CMB), and 
is often specified by giving the value of eg, which is the 
normalised mass variance in a sphere of radius 8/i _1 Mpc 
calculated from the PS when the model is extrapolated 
linearly to today (i.e. z = 0). 

The PS thus shows the HZ form at small k, reaches a 
maximum at fc t w 0.2(h. _1 Mpc) — 1 and then interpolates 
between approximate power-law behaviours P{k) ~ k n 
with n varying from n ~ — 1 to an asymptotic (large k) 
value of n = — 3 |38j |. In practice here we will not work, 
for our simulations of molecular dynamics, with the full 
PS described in Eq. l|57|l : our simulations are of a size 
which does not allow us to resolve the numerous different 
scales in this expression. We use instead a simplified 
version of this PS which retains its essential qualitative 
features: 



P(k) 



Nk 



1 + (Ak) a exp (k/k c 



(58) 



with the maximum k t chosen well inside the simulation 
box. 

Following the discussion in Sect. [H]we seek to produce 
a discrete distribution with PS S(k) given by Eq. 115(1 
with 



D(k) 



1 



-k 2 /2k 2 d 



(59) 



We note that, with this choice for the function D(k), 
the upper bound on kd, taking in Eci. (|18f) £ c (0 = and 
using the condition h(r) > —1, is: 



k d < \/2^(n ) 1/3 « 1.55/a, 



(60) 



where we have used the definition of a given in Eq. ((2311 . 
By increasing no sufficiently one can represent the con- 
tinuous model up to a desired k. 

In the first subsection below we will present an example 
of a HZ spectrum generated with a simple 1 /r 2 potential. 
In the following subsection we present the method using 
the simplified PS of Ea. 1)58(1 . while in the last subsection 
we give the potential which should allow the generation 
of the "realistic" cosmological PS of Eq. 1(57(1 . 



A. The HZ spectrum 

We consider just the "primordial" part of the PS with 
the HZ behaviour P(k) ~ k. As mentioned in Sect. IIIII 




md, r=o.6 
hnc, r=o.6 
md, r=i.6 
HNC, r=i.6 



ka 



FIG. 3: The PS of a 1/r 2 OCP for two different values of 
the coupling parameter F'. Excellent agreement is observed 
between the predictions from the HNC and MD in the range 
where they overlap. For the weak coupling case the HZ form 
for the PS S(k) oc k is clearly evident. The units are nor- 
malised to the ionic radius a. Note that the plot is on a 
linear-linear scale. 




md, r=0.6 
HNC, F=0.6 
MD, r = i.6 
HNC. r=i.6 



r/a 



FIG. 4: The correlation function in real space for the 1/r 2 
OCP for the same cases as in the previous figure. For 
the smaller coupling one has anti-correlation at all scales 
(g{r) < 1) while for the larger coupling one sees, just as 
in the standard OCP, the correlation appear with the first 
neighbour (which becomes more localised as the temperature 
is lowered). 



it has been shown in [Dj by a simple screening argument 
that one expects to obtain such a PS at small k in a 
modified OCP with interparticle potential v(r) = 1/r 2 . 
To verify this expectation we have used both the HNC 
and molecular dynamics as described above. In Fig. |3| 
the results for the PS are given for each case, and in Fig. 
0]the correlation function in real space. Because the po- 
tential is still a pure power-law the phase space is, as for 
the standard OCP, one dimensional and may be charac- 
terised by a single dimensionless parameter analogous to 
that for the OCP. We make the obvious generalisation of 
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the definition in Eq.J33J: 

T' = f3{Zef/a 2 . 



(61) 



The results from the HNC are valid in the infinit vol- 
ume limit and show very good agreement with the pre- 
diction of the asymptotic form for both the PS and the 
correlation function given in The range of these 

behaviours is, as expected, greater for smaller values of 
the coupling, and the linearity of the PS in particular is 
clearly visible in this case. We have checked also that one 
recovers the characteristic behaviour of the correlation 
function at large scales (g(r) — 1) ~ — 1/r 4 , which is also 
that of cosmological models with this PS (see 0,0])- 

The simulations of molecular dynamics were performed 
in the micro-canonical ensemble with the methods de- 
scribed in Sect. |IV] above, with N=4000 particles [39j . 
This corresponds to a simulation box with side of length 

L = (^§-N) 3 a as 25.6. Over this limited range very 
good agreement is seen with the results from the HNC 
in all cases, with some remaining statistical fluctuations. 
The units of time used in the simulations is r = 
with u>p — AimQ(Ze) 2 /m. To ensure good conservation 
of energy we have used a time increment of typically 
At ~ 10~ 2 t, which leads to fluctuations of ~ 1CP 7 in 
the energy. The system evolves for 10 5 r times steps, at 
which point it has reached thermal equilibrium. Then the 
PS and correlation functions are computed over many re- 
alizations of the system. By the ergodic principle this is 
equivalent to performing an ensemble average. Each re- 
alization is thus a configuration of the system at each 
time step. We compute the average in all the simula- 
tions over 50000 time steps, which leaves only very small 
fluctuations about the average. 



B. CDM-type spectra: simple model 

Let us now consider the spectrum (|58|) . We have seen 
that the small k part of the spectrum can indeed be pro- 
duced by a repulsive 1 /r 2 potential. We will now use the 
inversion of the HNC to determine the short range poten- 
tial which needs to be added to produce the modification 
at larger k described by Eq. (JSSJ. 

This input continuous model PS contains four free pa- 
rameters. We define S(k), the PS of the discretisation, 
using Eq. I|15|) with D(k) as in Eq. H59|) . and kd ex- 
pressed in terms of a by its bounding value in Eq. (|60|1 . 
In discretising we thus introduce just one independent 
length scale, fixed by the particle density no, and which 
we identify by convention with the ionic radius a de- 
fined in Eq. (|23[l . To fix the model being discretised 
we must specify S(k) by giving the three dimcnsionful 
quantities (Af, A and k c ) in units of a. The parame- 
ter a in Eq. (|58|l fixes the slope of the power spectrum 
P(k) ~ fc 1_a beyond the "turn-over" from the small k 
behaviour P(k) ~ k. As discussed above relevant values 
are in the range 2 < a < 4, and we will consider here the 
case a = 3 (i.e. P(k) ~ k~ 2 beyond the turn-over). 



The cut-off parameter k c is required simply to make 
the theoretical one-point variance (given by the integral 
of the PS) finite. We will take it, however, to be much 
larger than 1/a, so that it will play essentially no role in 
what follows other than to regularise Fourier transforms 
at small separations (e.g. in determining the theoretical 
correlation function £ c (r) associated to Eq. (J5SJ). The 
numerical results below used k c — 2.7/a. 

We thus have essentially just the two parameters A 
and Af to specify in terms of a. The former fixes the 
scale of the maximum (or "turn-over) in the PS: 



1 



1 



(a-l) 1 /^ A 



(62) 



When wc simulate the molecular dynamics we have a fi- 

1 /3 

nite box of side L = (-jW) a, and thus a finite range, 
from the inter-particle distance to the box size, in which 
we can attempt to represent the theoretical correlation 
properties. In k space this is the range kf < k < kd, 
where kf = 2-k / L is the fundamental mode and kd the k- 
space scale beyond which S{k) necessarily deviates from 
the input P(k). We thus choose A so that k t falls in this 
range. Once this choice is made Af then fixes the over- 
all normalisation. As discussed in Sect. Ill Dl the relevant 
dimensionless parameter to characterise the amplitude of 
the fluctuations is noP(kd). The natural choice, which 
we will use, is noP(kd) ~ 1 i.e. the theoretical power is 
roughly at the Poisson level at the inverse of the inter- 
particle distance. 

We will consider here a simulation of 8788 particles i.e. 
L w 33.3a and kf ~ 0.19/a. We choose kt = 0.52/a so 
that we have a small range of wavenumbers in which the 
PS is linear in k inside the box. For the chosen value of 
a = 3 this gives A 1.51a. Finally we take Af = 30a 4 
which corresponds to noP(kd) = 0.47. In Fig.[5]the pair 
P(k) and S(k) corresponding to these values are shown, 
and in Fig. [5] the corresponding correlation function pair 
£ c (r) and h(r). 

As discussed in Sect. Ill 51 the theoretical PS P(k) can 
be interpreted as that of a continuous mass distribution 
obtained by a physical smoothing with a function, whose 
Fourier transform is given by 



\w kd , no (k)\- 2 = i+ 



(l + (Afc)")(l-e- fc2 / 2fc ') 
noAfk 



exp(fc/fc c ). 

The corresponding real space smoothing function, calcu- 
lated numerically for our chosen S(k), is shown in Fig. 
[7| It decays at large separation faster than 1/r 4 , and is 
thus a localised smoothing in the sense we discussed in 
Sect. [H] It has, however, the rather unsatisfactory fea- 
ture of oscillating through negative values, albeit when 
the amplitude is already very small. We could, in princi- 
ple, remedy this by making a slightly different (but more 
complex) choice of D(k), but we do not anticipate that 
it should cause any significant change in our results. 

Given our chosen S(k) we can now use the inverted 
HNC equation Eq. I|39|l to determine the combination 
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— P(ka) 

— S(ka) 



2 

ka 



FIG. 5: The theoretical PS P(k) and its discretisation S(k), 
for the parameters choices given in the text. They differ as 
described by the chosen smoothing function D(k), for fca>l, 
with S(k) going asymptotically to 1/no, where no = 3/(4-7ra 3 ) 
is the particle density for the discretisation. 




FIG. 6: The theoretical input two-point correlation function 
£ c (r) and that of the corresponding discretisation h(r). 





FIG. 7: The window function in real space \Wk d ,n (r) 



FIG. 8: The correlation function, discretised correlation func- 
tion, and interaction potential obtained by the inversion of the 
HNC for a PS as given in Eq.^H), with JV = 20 a 4 , A ft! 0.69 a 
k c — 2.7 /a and kd = 1.55/a. 



j3v(r). The results are shown in Fig. |SJ For r/a > 5 the 
potential approaches the 1 /r 2 form corresponding to the 
small k behaviour of the PS. For r>4a the correlation 
function is negative (i.e. the system is anti-correlated) 
while it becomes positive for smaller scales. Correspond- 
ingly the potential becomes attractive at approximately 
this scale. At a smaller scale r < lsim2a, we see that the 
continuous correlation function £ c (r) (i.e. the Fourier 
transform of Eq. JHHJ) begins to deviate from the two- 
point correlation function h(r) of the discretisation. 

As anticipated, given the characteristics of S(k) which 
describes a positively correlated system at small scales 
(associated with the negative power law behaviour for 
k > k t ), the potential v(r) is attractive at small scales. 
Such a potential does not have a well defined thermal 
equilibrium [22^ . To define a system whose thermal equi- 
librium gives the desired S(k), we therefore add a re- 
pulsive core at some scale. We take a core of the form 
v c (r) = 0.2a 10 /r 12 (in units Ze = 1) which has a negligi- 
ble effect on the potential at and above the scale a. We 
check that this full potential, used in the direct HNC, 
indeed reproduces the input S(k) to a very good approx- 
imation. This is also a good consistency check on the use 
of the HNC for the chosen parameter values, in particu- 
lar the normalisation given by the chosen value of M. If 
J\f were too large we would no longer be in the regime of 
weak correlations for which the HNC equation is valid. 

We finally perform a simulation of molecular dynamics 
with the determined potential to obtain configurations of 
points with the PS desired. The HNC equation gives the 
potential times the temperature 0v(r), so we must choose 
a temperature. Any choice which keeps us in the regime 
of validity of the HNC should be appropriate. We make 
the simple choice j3 = a 2 (in units Ze — 1 i.e. in which 
the large separation limit of the potential is always simply 
v(r) = 1/r 2 ) which gives r' = 1 in Eq. (|61fl where, as we 
have discussed, the HNC gives very accurate results for 
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FIG. 9: The PS measured in a simulation of the molecular 
dynamics of 8788 particles for the potential shown in the pre- 
vious figure. Also shown is the input PS i.e. the PS of a 
system at equilibrium with this potential as calculated in the 
HNC. 
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FIG. 10: The real space correlation function for the same 
cases as in the previous figure. The abrupt break to anti- 
correlation at about r/a ~ 0.25 comes from the hard core 
introduced to ensure stability of the system. 

the 1/r 2 potential [iof . 

In the simulation we use this temperature to fix the ini- 
tial velocities, taking them to be random, and Gaussian 
distributed, with the corresponding variance. Typically 
the system thermalises at a slightly different tempera- 
ture (as there is also a change in potential energy relative 
to the initial configuration). By trial and error one can 
converge on the initial dispersion which gives the desired 
final temperature. For these simulations this involved 
typically a few trials. We will return to this point briefly 
in the final section below. 

In Fig. [5] and Fig. ^|are shown the results for molec- 
ular dynamics simulations with the potential determined 
from the HNC as described above (and for a final tem- 
perature s.t. (3 = 1/a 2 )- The agreement is satisfactory, 
in the sense that we reproduce quite well the correlation 
properties of the discretisation we have defined, in both 
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FIG. 11: Correlation functions and interaction potential ob- 
tained for the cosmological CDM spectrum described in the 
text. 



real and reciprocal space. The effect of the hard-core is 
visible in real space in the abrupt break in the correla- 
tion function, and it has a negligible effect above the scale 
r w 0.2a below which it dominates |4lj . 



C. CDM-type spectra: realistic model 

We consider finally the determination of the poten- 
tial which should reproduce the cosmological PS 1)57(1 . 
with the parameters of a currently standard cosmologi- 
cal simulation (e.g. like that taken as initial condition 
in the simulations of the VIRGO consortium 9] ) . To do 
so we must choose, in units of physical length, the scale 
a characterising the desired discretisation. For our ex- 
ample, we choose this scale by supposing we have the 
same physical density of point no as in some typical 
simulations of the VIRGO consortium: we suppose that 
we have the particle density corresponding to 256 3 par- 
ticles in a cubic box of side 239.5h~ 1 Mpc. The gives 
a as 0.58h _1 Mpc. We take our initial time to corre- 
spond to red-shift z = 50, and fix the normalization of 
the model at this time by the prescription that, using 
the extrapolation of linear theory, one obtains today (at 
z = 0) cr 8 = 1. This corresponds to a normalisation 
such that erg = l/(z+ 1) = 1/51, and the normalization 
factor is then Af = 29381 (h^Mpc) 4 . The discretisation 
scale kd introduced is chosen again at the bounding value 
kd ~ 1.55/a. In Fig. ^2 are shown the correlation func- 
tions and the potential for this case. We note the same 
1/r 2 behaviour at large scales, but the potential is more 
complicated at small scales due to the oscillations in the 
direct correlation function c(r) . By simulating the molec- 
ular dynamics with this potential with a sufficiently large 
number of particles, as we have done for a smaller num- 
ber of particles for the simpler cases, we should obtain a 
discretisation of the model with the properties Eqs. I|21|) 
and 
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VI. OTHER ISSUES AND CONCLUSIONS 

We have presented a new method to generate discrete 
distributions with desired two-point correlation proper- 
ties. It provides a promising alternative to the standard 
method for generating initial conditions for N-body simu- 
lation in cosmology used, which involves displacing parti- 
cles in a prescribed manner of a perfect lattice (or, some- 
times, "glassy" configuration). As discussed in detail in 
[Io| this latter method usually represents well the input 
theoretical PS in Fourier space at wave-numbers below 
the Nyquist frequency, but produces in real space a sys- 
tem with correlation properties which are a mixture of 
those of the initial unperturbed lattice configuration and 
those of the theoretical model. One obtains, in partic- 
ular, a two-point correlation function which is a rapidly 
oscillating function up to very large separations, which is 
a behaviour completely different to that of the theoretical 
model. In comparison the interest of this new method is 
thus that it can give (by construction) a faithful rep- 
resentation of the two-point statistical properties of a 
CDM-like spectrum in both real and Fourier space. In 
particular, in the examples we have considered, the cor- 
relation function in real space converges very well to the 
theoretical input correlation function at a scale slightly 
larger than a. Increasing the number of particles for the 
same physical size of the system, the interparticlc dis- 
tance a, and thus this scale, diminishes. In Fourier space 
the agreement is good (by construction) for wavenumbers 
less than roughly a factor of two smaller than the inverse 
of the scale a. 

In conclusion we now discuss the remaining steps which 
are required to apply the method fully in the context 
of cosmological simulations. Specifically we address the 
question of the generation of initial velocities required for 
such simulations, and then that of the modifications nec- 
essary to produce much larger distributions than those 
described here. Before doing so we first summarise briefly 
the method. 



A. Summary of the method 

The principal steps involved in going from an initial 
input theoretical PS P(k) to an output configuration of 
initial conditions are the following |42(: 

• 1. Determination of the function S(k), the PS of 
the discretisation of P{k) which we wish to gener- 
ate. As discussed in Sect. [H] this choice is neces- 
sary as the PS of the point distribution cannot be 
equal to P{k), because continuous and discrete PS 
have intrinsically different asymptotic properties at 
large k. The relation between the two is given by 
Eq. (|15|) . which, requires the choice of a smooth- 
ing function D{k). We have considered the simple 
choice D(k) = 1 — e~ k / 2fc d, and have taken kd as 
given by its upper bound as in Eq. (|60|l . The pas- 



sage from P(k) to its discretisation S(k) involves 
therefore the introduction of only one parameter, 
a length scale which we have chosen to take as 
a, determined by the particle density no through 
Eq. 

• 2. Choice of parameter values in S(k). Given the 
functional form of S(k) one must give the specific 
model by choosing the values of the parameters. 
This means simply that the length scale a must 
be associated to a physical scale of the model (or, 
equivalently, that the dimensionful parameters in 
P(k) be specified in units of a). In making this 
choice of parameters there are two essential points. 
Firstly, the system will be simulated with a finite 
number of particles N, and thus a finite box (of side 

L = (^-N) 3 a). The parameters choices should 
evidently mean that the important characteristic 
scales in P(k) fall well inside the range kf < k < kd 
(where kf is the fundamental mode of the box). 
Secondly the overall amplitude of P(k) should be 
specified by a choice of the dimensionless parameter 
noP(kd), and should make this parameter of order 
unity. 

• 3. Determination of the two-body potential. Given 
a specific S(k) the inversion procedure using the 
HNC, described in detail in Section UTTl is used to 
determine /3v(r) (where v(r) is the two-body inter- 
action potential) . If this associated potential does 
not give a stable thermal equilibrium, a hard core 
is added well below the scale a to ensure stability. 
Using this full potential (i.e. obtained by adding 
the core) one verifies that the direct HNC give the 
desired S(k). This is also a check on the applica- 
bility of the HNC at the chosen amplitude. 

• 4. Generation of realizations of the discretisation. 
The molecular dynamics of N particles interacting 
with the potential v(r) is simulated, starting from 
initial velocities corresponding to a chosen temper- 
ature. A simple appropriate choice of the latter is 
(3 = 1/a 2 (in units in which the long-range part of 
the potential is v(r) = 1/r 2 ). If the system ther- 
malises at a slightly different temperature, the ini- 
tial velocities can be appropriately modified by trial 
and error. 



B. Generation of initial velocities 

In the initial conditions of cosmological simulations of 
gravity, velocities must evidently also be ascribed to the 
particles of the distribution. We will now explain that 
these velocites can be generated for our point distribu- 
tions using the same method as in the standard algo- 
rithm. 

In cosmological models small initial density perturba- 
tions evolve, in linear theory, as a combination of a grow- 
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ing and a decaying mode. The choice of velocities which 
is almost always appropriate corresponds to the selection 
of only the growing mode. In practice this is implemented 
using the Zeldovich approximation 0, Q , which gives the 
motion of a fluid element in a self-gravitating fluid as 



x(q,t) = q + f(t)u(q) 



(64) 



where x is the (Eulerian) position of the fluid element, q 
its Lagrangian coordinate and f(t) is a function of (cos- 
mic) time t. In an expanding universe x is the comoving 
coordinate of the fluid element, related to its physical co- 
ordinate r by r = a(t)x. We will suppose that q is the 
(comoving) position of the fluid element at a time to i.e. 
/(*o) = 0. 

Using Eq. I|64[l in the continuity equation, expanded to 
linear order in the displacement field u, gives 



p(x,t) = p(q) l-/(i)V,-tf 



(65) 



where p(x, t) and p(q) are the mass densities in Eule- 
rian and Lagrangian coordinates respectively. The pecu- 
liar gravitational field, defined by g = f — -r obeys the 
Poisson-type equation 



V r -9 ■ 
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p(r,t)-- 3 



(66) 



sourced by the density fluctuations in matter about the 
mean density (po is the constant mean comoving mat- 
ter density of the universe, equal to the mean density in 
physical coordinates when a — 1), and is given, using 
Eq. lEH), by 



g = a(t)[f(t) + 2H(t)f(t)]u(q) 



(67) 



where H{t) = a/a. Substituting this relation in Eq. l|H6")l 
at t = to (when x = q), and then in Eq. (|65J) . it is easy to 
verify that, to linear order in the density perturbations 
about the mean comoving density po, 



p(x, t) - po 
Po 



^Gp f(t) 



f(t )+2H(t )f{t ) 



p(q) - Po 



(68) 

where we have chosen a(to) = 1. 

The evolution, on the other hand, of density fluctua- 
tions in Eulerian linear theory is given by 



p(x, t) - po = D , t \f>(Q) - Po 



Po 



Po 



(69) 



where D(t) is a solution of 



D(t) + 2HD(t) - ^£°D(t) = , (70) 
a 6 

normalised to D{to) = 1. Thus taking f(t) = D(t) — 
D(to) in Eq. (EHJ, one obtains Eq. ^ i.e. from Eq. ^ 
with f(t) — D{t) — D(to) one recovers Eulerian linear the- 
ory for the evolution of small amplitude density fluctua- 
tions. In particular choosing the growing mode solution 



(i.e. D + (t) = (i/io) 3 for a flat matter dominated cos- 
mology) gives fluctuations evolving in the pure growing 
mode as usually desired. 

Using Eq. Il64|) the physical velocity of a fluid element 
is f = Hf + v , where 



v = a{t)f(t)u(q) 



(71) 



is its peculiar velocity. Now through Eq. (|67|l it follows 
that 



/(*) 



f(t) + 2H(t)f(t) 



(72) 



To fix the initial velocities in our case we therefore sim- 
ply use the Zeldovich approximation as given by Eq. 16411 . 
the q being identified with the initial positions of our 
particles in the initial configuration corresponding to the 
required initial time to. Eq. I|72() evaluated at t = to and 
with f(t) = D + (t) specifies the velocities, in terms of the 
peculiar gravitational field g at each point in the initial 
configuration. 

It is perhaps useful to clarify the relation of this use of 
the Zeldovich approximation to the usual one in setting 
up initial conditions for cosmological NBS. To do so we 
write Eq. l|6^|l as 



x(R,t) = R + D + (t)u(R) 



(73) 



by a change of the Lagrangian variable R = q — u(q) 
[5|. One has D + (t — > 0) = so that the density in the 
coordinates R is that of the universe as t — * 0. In the 
standard method of setting up initial conditions in cos- 
mological simulations, these initial positions R are taken 
to be those of points in a "uniform" configuration: usu- 
ally a lattice, or sometimes also a "glassy" configuration 
0. The initial conditions (at t = t ) are then generated 
by moving the points to q = R + u(R) with the displace- 
ments u(R) given through Eq. Q67JI , the gravitational field 
being evaluated at the position R by a realization of the 
desired gravitational potential field (which in turn is de- 
termined by the theoretical power spectrum through the 
Poisson equation) . This is accurate to linear order in the 
displacement field. 

The crucial point which allows us to use the Zeldovich 
approximation is that this approximation does not re- 
quire, as is supposed in its usage just described in the 
standard algorithm for generating initial conditions, that 
the evolution of the fluid start from a perfectly homoge- 
neous configuration. Indeed the Zeldovich approxima- 
tion is simply an asymptotic solution of the evolution of 
a self-gravitating fluid away from a generic initial config- 
uration with some initial velocity and gravitational field 
[5J. More specifically all that is required to recover the 
linear amplification of density fluctuations which follow 
from the Zeldovich approximation, as described in the 
first part of this section, is that the density fluctuations 
are small. In this respect, we note that one can actually 
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consider the standard algorithm as a two step process in 
complete analogy to our method. In the first step parti- 
cles are displaced off a lattice (or other "uniform" parti- 
cle configuration) in order to produce a distribution with 
a power spectrum resembling as closely as possible that 
of the theoretical power spectrum P(k). This does not 
in fact require the use of the Zeldovich approximation, 
but only the continuity equation. Indeed from Eq. I|73l) 
with D(to) — 1, and assuming p{R) — po, the continuity 
equation gives 

Pm_£l = ^ R . a ( R) . (74) 

Po 

From this it follows that it is sufficient to take u(R) to 
be the gradient of a scalar field 4>(R) (as any rotational 
component for u will not generate density fluctuations) . 
Through Eq. we have that | </>(£) | 2 = P(k)/k 4 , and 
thus that 4>{R) is just proportional to the gravitational 
potential. The procedure is thus precisely the standard 
one. The second step then consists in generating the ve- 
locities, and one can proceed exactly as we have described 
for our method, using now explicitly the Zeldovich ap- 
proximation and, specifically, Eq. I|72[) with the gravita- 
tional field g given at each point in the perturbed con- 
figuration (i.e. at the Lagrangian coordinate q). To the 
linear order in the displacements this is identical to the 
standard prescription. 

C. Generation of larger distributions 

To test the proposed method here we have generated 
particle distributions with at most of order ten thousand 
particles, a number typical in statistical physics for the 
simulation of the molecular dynamics of systems of the 
type we have used (and for which the numerical codes 
have been conceived). It may be interesting and useful 
to evolve configurations of this size generated in this way, 
and by the standard method, under their self-gravity and 
to do a comparative study of this evolution. However, to 
use it as a truly alternative method in full cosmological 
simulations, it will be necessary to generate much larger 
configurations: in this context simulations typically now 
involve, at a minimum, hundreds of thousands of parti- 
cles. 

The Ewald summation method we have used here can 
be shown to scale (for a fixed precision on the cal- 
culation of the energy and forces) as TV 3 / 2 . Given that 
the largest simulations we have reported here, of approx- 
imately nine thousand particles, ran in about one day on 



a typical current PC, we could simulate only a slightly 
larger number in a reasonable time with our code. There 
are, however, well documented ways of speeding up the 
Ewald sum method by calculating the reciprocal space 
part of the sum on a grid. One can then employ rapid 
FFT algorithms, which scale as N log N |25|, with the 
direct space part of the sum scaling as N. There are nu- 
merous variants (e.g. P3M, PME, SPME), which differ 
primarily in how the charges (or masses) are assigned to 
the grid. Going beyond the Ewald method, or in com- 
bination with it, an even faster algorithm is the Fast 
Multipole Method (FMM) (see e.g. (26|). It is based 
on a decomposition of the system into sub- volumes, with 
the interaction between the charges contained in a given 
volume with those in a sufficiently distant sub-volume 
being calculated with a multipole expansion. Most of 
these techniques are of course known in the context of 
cosmological NBS (e.g. pj, 0) where they have been 
applied to gravity. Note, however, that their use if not 
limited to this case, and in particular does not require 
that the potential necessarily obey a local equation like 
the Poisson equation as for the 1 jr potential. 



D. Other remarks 

Another possible improvement concerns the choice of 
initial velocities in our simulations of the molecular dy- 
namics. When introducing the potential calculated with 
HNC we implictly choose the equilibrium temperature 
before performing the MD simulation. Thus, as ex- 
plained above, we put the initial velocities to get the 
chosen final temperature. The problem is that we do not 
know a priori how the system is going to reach equilib- 
rium and it is necessary to do trials with different initial 
velocities until the desired equilibrium temperature is at- 
tained. A solution to this problem would be to modify 
the MD to work in the canonical ensemble (in which the 
temperature is fixed) rather than in the micro-canonical 
ensemble as we have done here. 

We remark finally that we have not attempted here 
to address at all the very interesting underlying question 
motivating this work of the importance of the accuracy 
and nature of the discrete representations of cosmologi- 
cal IC in the dynamical evolution under gravity of these 
IC. Our purpose here has been to provide an alternative 
method for representing IC which should allow such ques- 
tions to be addressed with a new method, and hopefully 
answered, in due course. 

We thank J.M. Caillol, A. Gabrielli and F. Sylos Labini 
for useful discussions and comments. 
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